Columnar versus smectic order in systems of charged colloidal rods 
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We study the stability of inhomogeneous liquid crystalline states in systems of monodisperse, stiff, 
charged rods. By means of a bifurcation analysis applied to the Onsager free energy for charged rods 
in strongly nematic states, we investigate nematic-smectic and nematic-columnar instabilities as a 
function of the Debye screening length . While the nematic-smectic transition clearly pre-emts 
the nematic-columnar one in the regime of strong screening (i.e. small k,~^) a marked stability 
of hexagonal columnar order is observed at larger screening lengths. The theoretical results are 
substantiated by Brownian dynamics computer simulation results based on the Yukawa-site model. 
Our findings connect to experiments on tobacco mosaic virus rods in particular but might be relevant 
for soft rod-like mesogens in strong external directional fields in general. 
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I. INTRODUCTION 

Over the past decades, much research effort has been 
devoted to understanding the liquid crystal phase behav- 
ior of hard non-spherical colloidal particles, particularly 
in density functional theory and computer simulations. 
The theoretical approach to lyotropic liquid crystal for- 
mation has been initiated a long time ago by Onsager— 
with his classic paper on the isotropic-nematic transi- 
tion of infinitely thin rods. This theory shows that re- 
pulsive interactions alone can lead to long-range orien- 
tational (nematic) order, disproving the notion that at- 
tractive interactions are a prerequisite. The full phase 
diagram for hard spherocylinders with variable aspect ra- 
tio has only fairly recently been mapped out completely 
by means of extensive computer simulations^i^. A no- 
table feature is that nematic order is only stable for suffi- 
ciently large aspect ratios while isotropic systems of short 
rods tend to freeze directly into a crystalline state with 
three-dimensional long-range positional order. As to the 
nematic state, additional high-density phase transitions 
may occur involving partial freezing into liquid crystals 
with long-range spatial inhomogeneities in either one di- 
mension, giving a smectic A (SmA) phase, or two di- 
mensions leading to columnar (C) order. The latter also 
plays a dominant role in the high-density phase behavior 
of plate-like colloids where a stable hexagonal columnar 
was found in simulations^'^ and experimental studies of 
model clay suspensions^'^. 

Complicating issues, often important in interpreting 
experimental results, include the effect of rod flexibilitji^i^ 
and size polydispersity^ . Flexibility plays an important 
role in systems of e.g. stiff polymers and linear micelles 
and its generic effect is a significant depression of ne- 
matic order compared to rigid particles^*'. Systems of 
self-assembling worm-like micelles may, under certain cir- 
cumstances, not show nematic order at all. Instead a 
phase transition from an isotropic fluid to a hexagonal 
columnar phase takes plac o^^'^^ . The stability of the 
columnar state is also strongly stimulated by the inher- 



ent length polydispersity which may be of the annealed 
form, where the size distribution depends upon density 
such as in micellar systems, or of the quenched type like 
for many colloidal systems where the size distribution 
is fixed by the synthesis procedure. A simple packing 
argument suffices to understand that rods with variable 
length do not easily fit into layers pertaining to smectic 
order and therefore prefer columnar order. The crossover 
to columnar order has been observed explicitly in binary 
mixtures of hard rods with two different lengthsi^iii. 

Another important factor in the phase stability of 
rodlike mesogens is the infiuence of the soft repulsive 
interactions. In the experimental situation, interac- 
tions between colloidal model rods with chemically mod- 
ified surfaces can never be rendered truly hard^*^'^^ and 
many important biomacromolecular systems such as to- 
bacco mosaic virus (TMV) and fd virus rods^^, or min- 
eral systemsi^ like goethite^^ and vanadium pentoxide 
(V205)^ rods are stabilized by electrostatic particle re- 
pulsions. The influence of these interactions on the 
isotropic-nematic transition had already been addressed 
by Onsager in his original paper. Later on, the effect of 
electrostatic 'twist', which disfavors parallel rod conflg- 
urations and hence destabilizes nematic order, has been 
worked out in more detail in a study by Stroobants et 
al?\ 

Much more challenging however is the question how 
the electrostatic interactions affect the stability of the 
inhomogeneous liquid crystal phases, in particular the 
smectic phase, formed at high densities. Theoretical at- 
tempts in this direction have been undertaken in Ref. 
[2^ and m. Based on a model which combines the con- 
cept of an effective rod thickness (to account for the ex- 
tent of the electric double layer) with a simple cell de- 
scription to treat the inhomogeneous states, Kramer and 
Herzfeld^^i^i were able to construct the phase diagram of 
parallel, charged rods. Their main conclusion, a generic 
charge-induced stabilization of the inhomogeneous liquid 
crystal states also shows up in the study of Graf and 
Lowen^^ based on an elaborate density functional the- 
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ory. 



II. ONSAGER FUNCTIONAL 



In this paper we present a similar study starting from 
the Onsager theory for freely rotating rods. Rather than 
attempting to construct the full phase diagram we will 
merely focus on locating symmetry-breaking instabilities 
of the nematic state towards the inhomogeneous phases. 
The main candidates are the smectic and the hexagonal 
columnar phases. The study of these type of instabilities 
has been pioneered a few decades ago by Mulder^^'^*^ and 
the preferential stability of smectic over columnar order 
in systems of freely rotating hard rods, first recognized 
in simulations^-, could also be established in theory^®<. 

Here, we will elaborate on the Onsager-type ap- 
proaches and extend the calculations towards charged 
systems in two steps. First, systems of parallel charged 
rods will be considered by employing a straightforward 
site model to account for the electrostatic end cap effects. 
Next, the restriction of parallel confinement will be re- 
moved and the influence of rotations and electrostatic 
twist shall be explicitly taken into account by consider- 
ing the pair potential for infinitely stretched line charges. 
Although the second virial theory is not quantitatively 
valid for our case (the weight of parallel rod configura- 
tions which drive the translational symmetry-breaking 
instabilities necessitates inclusion of higher virial coeffi- 
cients into the free energy), we expect the theory to give 
a reliable qualitative sketch of the competitive stability 
of smectic and columnar order for soft rods. 

In contrast to the previous theoretical studies of Refs. 
[22I and m, our calculations reveal a distinct crossover 
from smectic to columnar order upon decreasing ionic 
strength. These results may be helpful in interpreting 
experimental results on TMV rods at low ionic strength. 
Most importantly, we show that columnar stability in 
rod systems need not be induced by length polydisper- 
sity but may be brought about by the soft electrostatic 
interactions only. Moreover, it is shown that the region of 
manifestation of columnar order can be broadened con- 
siderably if the isotropic phase is suppressed. This can be 
achieved by applying a strong external directional field. 
Preliminary Brownian dynamics simulations based on a 
Yukawa-site model are also carried out and the results 
point to a marked stability of columnar textures at low 
screening conditions for asymptotically aligned rods, in 
agreement with theory. 

This paper is organized as follows. In Sec. II the gen- 
eral Onsager functional is introduced and the bifurca- 
tion analysis is outlined in Sec III. The subsequent sec- 
tions deal with the explicit calculation of the appropriate 
Mayer kernels which form the necessary ingredients of 
the analysis. This is done first for parallel charged rods 
in Sec. IV and then for freely rotating hard rods (Sec. 
V) and charged ones (Sec. VI). Sec. VII is devoted to 
the Brownian dynamics simulations and all results will 
be combined and discussed in Sec. VIII. Finally, some 
concluding remarks will be formulated in Sec. IX. 



The starting point of our analysis is to construct a 
general free energy functional for an inhomogeneous sys- 
tem of N freely rotating rods with arbitrary mutual in- 
teractions in a volume V. The simplest approach is to 
take Onsager's well-known result^ for a virial expansion 
of the Helmholtz free energy F truncated after the first 
non-trivial term, which we may recast into the following 
functional 

PF[p] = J p(s){lnVp(s)-l}ds 

dspis) I ds'pis'Ms,s') (1) 

with f3~^ = ksT the thermal energy (fc^ is Boltzmann's 
constant and T temperature) and V the thermal vol- 
ume of a rod. The rod density distribution p is a func- 
tion of the generalized coordinates s — {r, fi}, indicat- 
ing position and solid angle, respectively and is normal- 
ized according to J dsp{s) = N/V. The first term in 
Eq. Il]) representing the ideal mixing free energy is exact 
while the second contribution contains the Mayer func- 
tion $(s,s') to account for the rod interactions on the 
approximate pairwise level. The Mayer function is re- 
lated to the direct rod pair interaction energy w via 



$(Ar; n, n') = exp[-f3w{Ar; Q, ft')] - 1 



(2) 



which depends on the centre-of-mass difference vector 
Ar = r' — r of a pair of rods and their orientations. For 
hard rods (3w is infinitely large if the particles overlap 
and zero otherwise. If the system is also homogeneous, 
the density distribution is independent of position and 
can be factorized into p ~ f{^l)N/V, where f{fl) is a 
distribution of orientations obeying J dflf{fl) — 1. Spa- 
tial integration of Eq. ([2]) then gives the relation 



J dAr<^{Ar;n,n') 



^GXCl 



(3) 



with Wexci = 2L^D\sm'y{n,n')\, the excluded volume 
of two thin hard rods with length L and diameter D 
(such that L/D ^ 1) at interrod angle 7. In gen- 
eral, Eq. ^ forms the basis of the classic Onsager-type 
theories for homogeneous systems of hard anisometric 
particles^S. For any given density N/V the thermody- 
namics of the system is fully contained in the orientation 
distribution function f{Q) (ODF). Its equilibrium form 
can be uniquely obtained for a given density by requiring 
the free energy to be minimal. Formally minimizing the 
free energy with respect to the ODF leads to the following 
general stationarity equation for homogeneous systems 



dfl'fin') J dAr^{Ar;n,n') 



/(17) = A/'-^exp 

^(4) 

with Af = J dri'exp[- • •] to ensure normalization. Upon 
increasing density, a change of shape of f{^l) from a 
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constant to a peaked distribution signifies an orienta- 
tional symmetry-breaking transition of the system from 
an isotropic state (with random rod orientations) towards 
e.g. a nematic one where the rods point in globally the 
same direction, defined as the nematic director. 

At higher densities, additional phase transitions from 
the nematic phase towards states with broken transla- 
tional symmetry such as smectic or columnar phases can 
be expected. A convenient way to approximately locate 
phase transitions of this kind is to apply a bifurcation 
analysis2^i2£ to the generalized Onsager free energy func- 
tional. 



III. BIFURCATION ANALYSIS 

In a smectic or columnar liquid crystal, the density is 
no longer constant throughout space but shows periodic 
spatial modulations in one dimension along the nematic 
director (smectic order) or in a two-dimensional Bravais 
lattice (hexagonal, cubic, et cetera) perpendicular to the 
director (columnar order). We may thus propose the fol- 
lowing Fourier expansion for the density distribution: 



(5) 



k l>0 



in terms of the orientation dependent amplitudes ai(f2) 
quantifying positional order along the various Fourier 
modes related to the wave number q — 27r/A, with 
A the typical spacing pertaining to the density mod- 
ulations. While for smectic order a single wavevector 
{k = 1) suffices, implementing columnar (or crystalline) 
order requires a superposition of different wave vectors 
{qijq2,-"} to reproduce the desired lattice. We will 
come back to this later. Identifying ao{^l) = pf{^) for 
the homogeneous system, we may simplify the above ex- 
pansion as follows (omitting the ^-summation for clar- 
ity): 



p{v,n)^pf{n) 



1 + ^ ai cos(Zq • r) 



(6) 



where we used that q = — q. The ODF f{Vl) is now 
taken to be that of the nematic reference state for all 
modes I. This means that any coupling between fluc- 
tuations in the spatial density and the orientations is 
neglected. A fully consistent calculation for the nematic- 
smectic bifurcation"^^ shows that the position-orientation 
coupling is marginally small for strongly aligned, slender 
rods on which we focus in this study. 

Inserting Eq. ([5]), truncated after Z = 1, into the free 
energy and expanding the free energy of the new inho- 
mogeneous (I) state with respect to the homogeneous ne- 
matic (N) one 5F = Fj — Fn up to quadratic order gives 
the free energy curvature 

s^(3F = alp^ [ dnf{n) [ dn'fin') 



The nematic state becomes locally unstable if the cur- 
vature vanishes, S'^f3F = 0. The above expression then 
simplifies to the bifurcation condition 



(8) 



where the hat denotes a cosine transform of the Mayer 
function: 

l>(q; 17, 17') / dAr$(Ar; 17, 17') cos(q • Ar) (9) 



which is the key input of the analysis. The bifurcation 
density is defined as the smallest non-trivial physical so- 
lution p* of Eq. ([5]) with associated wave vector q* . Once 
the instability has been located, additional information 
about the thermodynamic stability and the order of the 
phase transition can be inferred from a parametric ex- 
pansion about the bifurcation point, as outlined in Ref. 
|26| . In Appendix B we shall reproduce the analysis and 
derive a general expression for the Landau free energy 
of the inhomogeneous state with respect to the homo- 
geneous nematic system. This result will then be used 
to verify the thermodynamic stability of the smectic and 
columnar states. In the next sections we will derive ex- 
pressions for <I>(q) for the nematic phase focussing first 
on parallel charged rods. The theory is then extended to 
include the effect of rotations. 



IV. PARALLEL CHARGED RODS: 
HARD-CORE YUKAWA SITE MODEL 

To calculate the Mayer kernel for parallel charged rods 
we have to find a suitable route to include the infiuence 
of electrostatic end effects. These end effects are intri- 
cately difficult to quantify"^^ but are nevertheless essen- 
tial in our description since the formation of stable smec- 
tic/columnar liquid crystalline structures at high densi- 
ties is driven primarily by correlations between the end 
cap of one rod and the main (cylindrical) manifold of the 
other as embodied in the 0{LD^) contributions to the 
Mayer kernel. To make headway we consider a model 
in which a rod with finite length L is composed of an 
array of n 1 spherical beads with diameter D placed 
at equidistant intervals. Each bead i {1 < i < n) from 
one rod interacts with bead j from the other via a hard- 
core Yukawa (HCY) potential. The total rod potential 
depends only on the centre-of-mass distance s (expressed 
in units D) and is given by 



CXp[— K£>(Sij — 1)] 



if all Sij > 1 



if any Sij < 1 



(10) 
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with K the Debye screening constant and e a dimension- 
less contact potential e = [Z /nf{\B / D) / {l + nD /2f de- 
pending on the total rod charge Z (in units of the elemen- 
tary charge e) and the Bjerrum length \b = Pe^ /e^Sr 
which is determined by the temperature and the dielec- 
tric constant £,■ of the solvent. The quantity Sij represent 
the distance (in units D) between site i from one rod and 
i from the other. It is convenient to introduce cylindri- 
cal coordinates, so that s — {s^ cosw, s±^ smuj, S||}, with 
expressed in units D and S|| in units L. The site-site 
distance Sy is then given by 

Sy = y4TcZV^^)2[^Td(j^^, z,j-l,2,---n 

(11) 

Here, d = {L/D)/{\/n + \\/n — 1) is the dimensionless 
spacing between two neighboring sites on a rod, chosen 
such that the quadrupolar moment)^ of the rod matches 



where the spatial integration has been performed over the 
region outside the cylindrical manifold determined by the 
hard-core excluded volume of two parallel cylinders. The 
complementary integral over the inner region where the 
cylindrical cores overlap yields the hard-core contribution 
(denoted by "H" ) 

^h{Qc,Qs) = 27rLi?2 [ 2ds\\ f s^ds^'PncY 
Jo Jo 

xJo{QcS±_)co8Qss\\ 

= -AttLD'MQs)^^^ (15) 

with Jo (a;) = sinx/x a spherical Bessel function and Jj. 
a standard one. In addition we used that ^hcy = ^1 
for overlapping hard bodies. The result is identical to 
that of Ref. [2y. We remark that the excluded volume 



that of a homogeneous line charge of length L. 

The total transformed Mayer function for the HCY-site 
model is then given by: 

^HCY (q) = LD^ J ds^HCY (s) cos(I?q • s) (12) 

The generalized wave vector is given by q — {qcO^Qs} 
in terms of smectic (S) and columnar (C) components. 
Hence, we can write Dq ■ s = Qcs± cosw + Qss\\, with 
the dimensionless wave numbers 

Qs = 2. (A) , . 2. (^) (13) 

related to the lattice spacing A. 

Integrating over the angle lu then yields for the electro- 
static contribution (denoted by "Y" ) to the Mayer kernel: 



(14) 



between two spherical end caps ~ 0{D)^ is not resolved 
exactly by the integrations of Eqs. and (|15p and the 
Yukawa contribution in Eq. (|14p is therefore expected 
to be reliable up to 0{LD^), which suffices for slender 
rods. The integrations in Eq. can be carried out 

numerically by imposing appropriate cutoff distances for 
sj_ and S||. The bifurcation condition Eq. ([5]) is obtained 

by substituting l> + $y and f{Q) = 6{n). 

V. FREELY ROTATING HARD RODS 

Before embarking on our calculation of $ for freely 
rotating charged rods, we will first discuss the effect of 
rotations on the translational symmetry-breaking bifur- 
cations in systems of hard rods. A closed expression for 
$(q) for freely rotating, hard spherocylinders has been 
derived by van Roij^: 



/•OO pOG 

^y{Qc,Qs) = 2ttLD^ J 2ds|| J s±ds±<S>HCY{s±, s\\)MQcS±) cosiQssii) 

/OO pi 
2rfs|| J s^ds^^HCYisj_, s\\)JoiQcSj_) cos{Qss\\) 



r 



*ff(q;^^i,^^2) = -2i D\sinj\jo{Dq--v)jo\ -q-<vj_] jo[ -q-W2 



L 



L 



-LD'^jo [ -^q- wi 
LD^jo ( ^q- W2 



27rcos ( ^q- W2^ + 2 sin (^^q • W2 ) 1^(^^1,01) 

27r cos ( ^q ■ wi^ ^^i^^ - 2 sin (^^q ■ wi ) W{M2,^2) 
I 



(16) 



which contains all contributions up to 0{LD^). Cor- expressions are available, can be safely neglected for suf- 
rection terms of 0{D^), for which no closed analytical 
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ficiently anisometric rods L/ D ^ 1. The inner products 
contain the rod orientation unit vectors Wi and W2 with 
associated unit vector v = wi x w2/sin7. The function 
W is defined as 



tt/2 



t/2 



daji{M cos(a - cj))) (17) 



W{M, (j)) = 

with ji [x) = smx — x~^ cosx a spherical Bessel func- 
tion. 

Defining a third unit basis vector Ui = Wj x v the 
quantity Mi and the angle (pi are given by 



M, = i?9V(q-v)' + (q-a.)^ 
q Uj 



i = 1,2 



arccos 



^(q-v)2 + (q.u,)2 



(18) 



The result for perfectly aligned rods is easily recov- 
ered by replacing the orthonormal set of basis vectors 
{ui,v,Wi} by the set of Cartesian unit vectors. We 
may thus identify ui = U2 = {1,0,0}, v = {0,1,0} 
and wi = W2 = {0,0,1}. With help of the identity 
joix) — cos{x/2)jo{x/2) Eq. ([16]) then immediately gives 
back Eq. P^ . 

Let us now perform an asymptotic analysis of the 
Mayer kernel, valid for strongly nematic systems where 
the average deviation of the rod vectors from the nematic 
director is small. 



A. Smectic symmetry 

For the smectic case, density modulations occur along 
the nematic director fi fixed at n = {0, 0, 1} and 
hence q = (/^{0, 0,1}. The orientation of rod i on 
the unit sphere can be parametrized in terms of a po- 
lar angle 9i and an azimuthal one Lpi, so that = 
{sin 9i sin (pi , sin Oi cos Lpi , cos Oi } . Expanding all relevant 
inner products for small polar angles 9i up to leading 
order gives: 



^qw. 



Dq - V 



7lC2 



yef+^f^^2M2«)sA^ 



(19) 



with A93 ^ ^1 — ^P\- For large aspect ratios {L/D ^ 1) 
all contributions of O [{D /L)9i] become marginally small 
so that the following limiting values can be deduced : 
M, and (ZJq • v) ^ 0. Consequently, Ji{M)/M 
1/2 and W{M,(j)) 0. Using this in Eq. ^ together 
with the identity jo{x) — cos(a;/2)jo(a;/2) gives the fol- 
lowing asymptotic expression for the transformed Mayer 
function for the smectic symmetry: 

$H(q;^^i,f^2) - -21"" D\ sin jiQi, n^MiQ s /2) - 

27:LD^joiQs) (20) 



The latter contribution is simply the result for parallel 
rods, given by Eq. (fT5|) taking the limit Qc 0, while 
the first is the leading order correction term arising from 
the rotations. 

The bifurcation condition is obtained by performing 
an orientational average according to Eq. ([8]). Following 
Odijk^, we introduce the Gaussian trial ODE 



a 



exp [— ^aS^] 



if < 6* < f 



fcid) ^ J- < (21) 
' exp [-ia(7r - 9))^] if f < < tt 

which depends only on the polar angle in case of uniaxial 
nematic order. The variational parameter a is chosen 
such as to minimize the nematic free energy. It must be 
much larger than unity because the normalization factor 
is only valid in this limit. After minimization a attains 
a quadratic density dependence, a = Ac^/ir where c is a 
dimensionless rod concentration c = {tt / A) pL"^ D related 
to the volume fraction 77 via c = rjL/D. We may now use 
the following asymptotic result'^: 



((|sin7l))/o-(^ 



1/2 



(22) 



where the brackets denote an orientational average. With 
this, the bifurcation equation Eq. ([8]) finally becomes 



4j2(Qs/2) - SvjoiQs 



1 



(23) 



Note that the rotational correction term is independent 
of the density. 



B. Columnar symmetry 

Let us now turn to the columnar state where a two- 
dimensional modulatory density pattern _L fi is expected. 
In case of a hexagonal Bravais lattice, a combination of 
the following three unit wavevectors: 



/V3/2\ 

qi = I M ' q2 = I 1/2 1 



^3 = 




(24) 

correctly reproduces the six-fold symmetry of the colum- 
nar density modulations in the xy plane. We may now 
repeat the asymptotic analysis to obtain the leading or- 
der expressions for the inner products. For the qi-mode 
we obtain 



2qi-w. 

Dc[i ■ V 

M, 



^Qc^O^ sin (pi 



Q 



— 01 cos V?! -I- 6*2 cos (^2 \ 

y^9l+9l - 29i92 cosAyj / 
c (25) 



The expressions for the other columnar modes show ex- 
actly the same scaling with respect to 9 and L/D and 



6 



are not shown. A quick inspection reveals that the 
asymptotic analysis does not bring much relief since i) 
(q-v) ^ C(l) and ii) the combination {L/D)9i appearing 
in the first line also remains finite after having performed 
the oricntational average. It can be estimated with use 
of the Gaussian ODF Eq. pT|) and the resulting scaling 
relation reads {{L/D)\e\)f^ ^ r]'^ ^ 0{1). 

Unlike for the smectic case, the complicated nature 
of the orientation dependent terms precludes an analyti- 
cal calculation of their Gaussian averages. Nevertheless, 
the integration over the angles 9i and ipi can be carried 
out numerically without difficulty using standard Simp- 
son quadratureSS.. The bifurcation condition Eq. ([8]) for 
the columnar symmetry is given by a linear superposition 
of kernels 



k=l 



l,(/7i;fc'2,</52, 



/g 



(26) 



substituting the asymptotic forms Eq. 
and using the Gaussian ODF Eq. ([2T|) . 



3 in Eq. (H 



C. Bifurcation points 

The solutions of Eqs. (l23|) and ([26|) have been col- 
lected in Table I. To obtain more realistic volume frac- 
tions a density rescaling according to Parsons-Lee (PL) 
method^i^i^ has been carried out. The approach, not 
treated here, can be implemented quite simply by rescal- 
ing the density according to p — > pgcs{ii)^ and equiva- 
lently 77 -> 775cs(?7): where gcs(??) = (1- (3/4)?7)/(l-?7)^ 
originates from the Carnahan-Starling equation of state 
for a hard sphere fluid. All volume fractions mentioned 
throughout the rest of the paper are obtained from this 
treatment. 

It is clear that both N-S and N-C are destabilized due 
to the effect of rotations. The decrease of the smec- 
tic layer spacing can be intuitively understood from the 
fact that the rod length projected onto the nematic di- 
rector (the 'entropic' length) is smaller than the bare 
rod length due to the oricntational fluctuations. Simi- 
larly, the increased columnar spacing can be attributed 
to an entropic rod diameter larger then the bare one. 
The N-S pre-empting the N-C one is in accordance with 
simulations^"''*'' and experimental observations^-. In ad- 
dition, the location of the nematic-smectic transition pre- 
dicted here compares very well with the value « 0.418 
reported from simulations^. 



VI. FREELY ROTATING CHARGED RODS 

Looking back at the expressions derived for hard rods, 
in particular Eq. ()20|) . we may conclude that the asymp- 
totic Mayer kernel consist of two parts. First, a 'ref- 
erence' part for (near)-parallel configurations associated 





'7s 


Qs 


■nh 


Qh 


parallel rods 


0.575 


4.493 


0.945 


5.136 


parallel rods (PL) 


0.338 


4.493 


0.441 


5.136 


rotating rods 


0.792 


4.858 


1.540 


4.424 


rotating rods (PL) 


0.404 


4.858 


0.543 


4.424 



TABLE I: Overview of the nematic-smectic and nematic- 
columnar bifurcations for hard rods. (PL) refers to Parsons- 
Lee theory. 



with the O(LD^) contributions in Eq. (HH)) and, sec- 
ond, a part of 0{L^D) which constitutes a non- vanishing 
contribution in the limit of asymptotically strong align- 
ment. In this section we will first derive the leading or- 
der 0{L^D) correction due to electrostatic effects. These 
correction terms will then, together with previous results, 
be compiled into expressions for the total Mayer kernel 
for charged rods in near-parallel configurations. 

Following Ref. [2l| we start with introducing the fol- 
lowing form for the electrostatic interaction between two 
infinitely long charged rods at shortest distance x and 
mutual angle 7: 



Pweiix; Oi, ri2) = A 



e xp[-/.(|x| - D)] 
|sin7(rii,r22)| 



(27) 



If the rod charge is not too high, the prefactor A can 
be expressed in closed form within the Debye-Hiickel ap- 
proximation 



A 



{KDf{\BlD)Kl{nDl2) 



(28) 



with Ki{x) a modified Bessel function. This expression 
stems from the linearized Poisson-Boltzmann equation 
for an infinitely long, charged cylinder with diameter D. 
The key quantities here are the linear charge density \bi^ 
expressed in terms of unit charges per Bjrerrum length 
and the ratio of the electric double layer thickness and 
the rod diameter 1/kD. The denominator in Eq. ([27]) 
embodies the electrostatic 'twisting' effect whereby a rod 
pair is energetically stimulated to adopt a perpendicular 
configuration. Clearly, for strictly parallel configurations, 
the expression is ill-defined which means that a N-I bifur- 
cation analysis cannot be based solely on this potential. 

The electrostatic line charge contribution of the cosine 
transform of the Mayer kernel can now be written as 

$ei(q;f^i,!^2) = J d|x|{exp[-/3u;ei(|x|;Oi,f]2)]-l} 

|x|>D 

X cos(q • x) (29) 

If we neglect end effects, the shortest distance unit vector 
x is equal to v and the following parametrization x = 
Dtv, with 1 < t < 00, can be applied. Rewriting the 



7 



integral gives 

/oo 
dt{exp [-(A7|sin7|)e-«^*] -1} 

X cos{tDc[ ■ v) (30) 

with A' = Aexp[KD]. To solve the integral it is ad- 
vantageous to make the following substitution u(t) = 
{A' /\ sin j\) exp[—KDt] and to recast it into a complete 
Fourier integral (denoted by the tilde). This gives: 



duu ^{e " — 1) exp [— i(7ln(u/uo)] (31) 



with Um = u{m) and q = k ^q- v. The complex solution 
of the integral reads 



iq 
-Mo 



+ r(-ig, ui) - T{-iq) 



(32) 



in terms of the complete and incomplete gamma func- 
tions, r(z) and r(a, z) respectively. Because of the co- 
sine transform in Eq. ([30]) we only need the real part 
Re(K$e/). In the asymptotic limit, ui ~ A/|sin7| ^ 1 
and the incomplete gamma function becomes negligibly 
small so that it can be omitted in the remainder of the 
analysis. Writing out the complex functions one arrives 
at 



sin qkD 

q 



Re 



'r(zg) 



(33) 



To further approximate the latter term let us first note 
that for the smectic and columnar modes we have 



qs = (kL) iQs(q- v) < 1 
qc - {kD)-'Qc{ci-v)^0{1) 



(34) 



recalling that kD ^ 0{1) and L/D 1. For the smectic 
symmetry, we may approximate 



iq \ 



- T{iq) w 7B 
with Euler's constant to obtain: 



9< 1 



(35) 



= DjoinDqs) 



sin (^KDqs) 
QS 



-fE COS UDqs) \ (36) 



with D ^ D {l + {KDy^{\nA~ ln| sin7|)}. It is easy 
to verify that the electrostatic contribution vanishes in 
the hard rod limit {kD — > oo) as it should. Taking the 
limit ^ Eq. ([55)1 simplifies to: 



(In A — In I sin7| -|- ^e) 



(37) 



We may now compose the total Mayer kernel by replac- 
ing the thickness contribution jQ{Dq ■ v) in Eq. by 



Eq. ((37)) . The contributions depending on L are left un- 
touched since the line charge model, by definition, can 
only affect the interaction thickness of the rod. These 
considerations lead to the following form of the Mayer 
kernel of two line charges: 

$(q;r!i,l)2) = ~2L2i?eff(r!i, 1^2)1 sin 7|j2(gs/2) (38) 

valid for the smectic symmetry. It is similar to the first 
term in Eq. (j20p but with D replaced by an orientation- 
dependent effective thickness Dcs'- 

D^s{ni,rt2) = |l + ^ (In^ - In I sin7| + 7^) 

(39) 

which is exactly the same as the effective thickness show- 
ing up in the second virial coefficient of two charged rods 
in Ref. M- 

The next step is to perform the orientational average 
using Eq. (j2ip . The variational parameter a now follows 
from minimizing the nematic free energy for charged rods. 
The associated minimum condition reads^ 



a - /Cica^/^ /C2cai/2(2 - In a) = 



(40) 



with constants /Ci = tt^^/^ p + (kD)-^{2 In A + 37_e - 2)] 
and IC2 = n~^/^{KD)^^ . This equation has to be solved 
numerically for any given concentration c and electro- 
static parameters A and kD. One may verify the hard 
rod limit kD — > 00 to obtain an = 4c^/7r. The orienta- 
tion average of Eq. (|38p can be worked with the aid of 
Eq. and the following asymptotic result^i 

(( |sin7|ln|sin7|))/^ ~ (^y^' (ina - 2 -|- 7^5) (41) 

It is expedient to normalize the variational parameter 
a = a/ an, so that a approaches unity in the hard rod 
limit. Using this, the nematic-smectic bifurcation condi- 
tion Eq. ([8]) for charged rods finally becomes 

1 = -4j2(Q5/2)ci-i/2 

~8r]jo{Qs) + P^Y{Qs,0) 



(42) 



with given by Eq. (|M|) . Comparing this with the 
hard rod result Eq. (|23|) we see that the first contribu- 
tion is now implicitly dependent on concentration, albeit 
weakly, due to the nonlinear character of Eq. (|^(])) . 

For the columnar symmetry we have to retain the com- 
plex gamma functions in Eq. (|33p and the resulting ex- 
pression reads: 



4e/ = DjoiKDqc) 



' {cos (nDqc) Re [r (iqc)] 



■ sin (nDqc) Im [F (iqc)]} 



(43) 



With the aid of Eq. (|35p , it is easily shown that the con- 
tribution vanishes in the limit kD 00 (corresponding 
to qc 0). 
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We may combine this with the remaining contributions 
that depend only on the rod length L. The Mayer ker- 
nel of 0{L^D) for two line charges within the columnar 
symmetry can then be written in the following form: 

$ = -2L2i:>|sin7|jo (^Y^- wi^ jo (^-^q- 

xG(q;r!i,02) (44) 

where G now replaces the Bessel function Jq^Dci ■ v) in 
Eq. (HID 

G(q;r!i,l]2) - -iKDr^{cos{KDqc)Re[T{iqc)] 

+ sin {KDqc)lm[T{iqc)]} (45) 

and the reader may verify that this contribution is recov- 
ered from G in the hard rod limit, as it should. With Eq. 

the total Mayer kernel for the columnar symmetry 
is proposed to be of the following form: 

$ = l'L2^(q; ni,n2) + $LD2(q; ^^1,^2) + $y (o, Qc) 

(46) 

with $L2£i given by Eq. and $^£(2 the hard-core 

contribution of 0{LD^) from Eq. pB]) . both depend- 
ing intricately on the rod orientations. Since the latter 
contribution no longer pertains to the excluded volume 
of parallel rods, as was the case for the smectic symme- 
try, the leading order term from the Yukawa site model 
$Y is expected to be orientation dependent too. How- 
ever, getting access to this contribution requires a full 
numerical integration over all spatial and orientational 
variables of the HCY site model for freely rotating rod 
pairs. Resolving the Mayer kernel would then give a 9 di- 
mensional integration (including the site-site summation) 
which clearly is a formidable numerical task. Therefore, 
we shall rely on the parallel contribution. The approx- 
imation can be justified in part from the fact that the 
variational parameter a and hence the degree of nematic 
order becomes rather large for sufficiently large Debye 
lengths, as we will see later on. Finally, the nematic- 
columnar bifurcations can be computed by inserting the 
full Mayer kernel into the bifurcation condition Eq. (1^51) 
and performing a numerical averaging over the orienta- 
tional degrees of freedom using the Gaussian ODF. 

VII. BROWNIAN DYNAMICS SIMULATIONS 

To supplement the theory, Brownian dynamics (BD) 
simulations have been carried out for a system of point 
Yukawa-site rods. The corresponding interrod potential 
is virtually identical to Eq. pU|) with omission of the 
hardcore contribution for Sy < 1 (achieved by taking the 
limit D ~^ 0). This means that overlaps of the inner 
cylindrical cores are in principle allowed in our simula- 
tions. However these configurations are rare because of 
the significant energy penalty involved. Another differ- 
ence with the description in Sec. IV is that the rods are 



no longer fixed in parallel configurations but are allowed 
to rotate freely. Apart from the rod length and the Bjer- 
rum length, the relevant lengthscale of the point Yukawa 
model is the Debye length. The ratio of the latter two 
can be expressed as 

kXb = ^i7TZp{\B/Lf (47) 

in terms of the dimcnsionless rod concentration p = 
NL^/V. 

The simulations comprise a finite difference integra- 
tion of the Langevin equations for interacting Brownian 
macro-ions according to the scheme of Ermak^^. For a 
detailed exposition of the update equations for Yukawa- 
site rods the reader is referred to the paper of Kirchhoff 
et alM: The short-time self-diffusive behavior of the rods 
is characterized by two translational diffusion coefficients 
(one parallel and perpendicular to the rod axis) and a ro- 
tational one. All of these depend on the hydrodynamic 
aspect ratio of the rods^^ for which we take a value of 16, 
comparable to that of TMV rods-ii. 

Systems consisting of = 500 rods each with n = 13 
sites were simulated in a cubic box with periodic bound- 
ary conditions. The dimensionless concentration was 
fixed at a low value p — 4.0. Starting configurations 
were generated by randomly inserting rods in a parallel, 
non-overlapping configuration. Equilibration is then car- 
ried out until the system has reached a stable isotropic 
configuration, characterized by a vanishing nematic or- 
der parameter, defined as the maximum eigenvalue of 
the tensor 

Q ^ ^ E ^ ® - (48) 

i 

where (S) denotes a dyadic product, I the second-rank unit 
tensor and the brackets a canonical average. The next 
step in our simulations is to nematize the system by ap- 
plying an external directional field along the z-direction 
of the simulation box. This gives rise to an external po- 
tential energy per rod, given by 

= -ecos2(w, -z) (49) 

in terms of a dimensionless field strength ^ > (in units 
fc^r). Obviously, a stable nematic state could also have 
been obtained by increasing the density of the system. 
However, we found that the simulations become increas- 
ingly cumbersome for large densities due to the slow dy- 
namics and the formation of metastable, transient states. 
Moreover, at high packing fractions the results are ex- 
pected to be sensitive to the details of the model (in 
particular the number of sites per rod) and many time- 
expensive test runs need to be carried out. 

Once the directional field is switched on, the (instan- 
taneous) nematic order parameter is found to increase 
rapidly until a plateau value is reached after some time 
interval. The associated average nematic order parame- 
ter is found to be close to unity for ^ = 50 and ^ = 25 
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FIG. 1: (a) Bifurcation diagram for parallel charged rods. Plotted is the normalized Debye length kT^ /D versus the hardcore 
volume fraction 77. At {kD)~^ > 0.444 a stable nematic-columnar bifurcation pre-empts the nematic-smectic one. The volume 
fractions are obtained from Parsons rescaling. (b) Dimensionless wave number Q [Eq. (|13[l ] corresponding to the bifurcation 
lines shown in (a). 




FIG. 2: Simulation snapshots showing the projections of the 
rod centres-of mass (indicated by dots) onto the a;i/-plane per- 
pendicular to the field direction. Results correspond to var- 
ious values for the external field strength ^ [Eq. (|49p]. (a) 
parallel rods 00), (b) freely rotating rods at ^ = 50, (c) 

same for ^ = 25 and (d) ^ = 10. 



duce the intralayer pair correlation function, defined as*^ 

= ]^ ^EE '^(^ - 1^'^ X ^D© (V2 - r., • z)^ 

^ (50) 
with 5 the Dirac delta function and 8 the Hcaviside step 
function. It represents the probability to find a parti- 
cle at distance r from a reference particle within a slab 
of thickness L perpendicular to the director. Basically, 
gi provides information about positional order perpen- 
dicular to the nematic director and a profoundly peaked 
function is to be expected in case of a columnar or crystal 
structure. Second, the parallel pair correlation function, 
given by 




monitors liquid structure parallel to the nematic direc- 
tor and a peaked gp will appear in case of smectic or 
crystalline order. To make sure the results do not suffer 
from finite size effects, additional simulations have been 
carried out using N = 900 rods with n = 13 and no 
qualitative differences were observed. Furthermore, the 
pressure tensor^^ has been monitored in all cases and no 
evidence for spurious, non-isotropic stresses induced by 
the cubic periodic boundary conditions was encountered. 



VIII. RESULTS AND DISCUSSION 



(near parallel systems) and around 0.8 for the lowest 
nonzero applied field strength ^ = 10. The nematic state 
is then equilibrated further to allow for possible transla- 
tional freezing transitions. To detect smectic order (along 
z) or columnar order (perpendicular to z) we first intro- 



Before discussing the results we first have to specify the 
electrostatic parameters of the rod solution. For these 
we take typical values for TMV rods suspended in wa- 
ter (e^ = 78) at neutral pH and room temperature. The 
rod dimensions are D = 18 nm and L = 300 nm and 
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FIG. 3: (a) The intralayer and (b) the parallel pair correlation functions for various field strengths ^. 



the Bjerrrum length is As = 0.716 nm. The total rod 
charge is fixed at Z = 390 and the corresponding lin- 
ear charge density is set at As = 1 (i.e. one unit charge 
per Bjerrum length) which are reasonable values if the ef- 
fect of counterion condensation is taken into accoun^iSiii. 
These parameters justify the use of the Debye-Hiickel ap- 
proximation to determine the potential amplitude in Eq. 

In Fig. 1 results are presented for parallel charged rods 
with n — 30 HCY sites. Increasing the number of sites 
per rod did not lead to significant changes in the bifur- 
cation diagram. At high ionic strengths [{kD)^^ < 0.1] 
the rod charges are highly screened and the phase be- 
havior is virtually unaffected by the charge. A sta- 
ble nematic-smectic bifurcation pre-empts the nematic- 
columnar one. At lower ionic strength the transition 
shifts to lower volume fraction due to fact that the en- 
hanced pair interactions lead to an effective 'interaction' 
volume fraction higher than the bare one. A qualitative 
change of scenario occurs at even lower ionic strengths 
[{kD)~^ > 0.44] where a nematic-columnar instability 
pre-empts the nematic-smectic one and a stable colum- 
nar phase can be expected. In this regime, the extension 
of the double layers is significant and very little pack- 
ing is needed for the system to order. The generic sta- 
bilisation of positional order is in accordance with the 
results of Ref. [12 and [H. The enhanced interaction 
range also shows up in the bifurcation wave numbers in 
Fig. lb where the sharp decrease implies a columnar 
spacing spanning multiple rod diameters. The jump at 
{kD)~^ « 0.28 is a consequence of the fact that there 
are two independent bifurcating solutions (i.e. a high-Q 
branch and a low-Q one) whose densities merge at that 
point. 

The stability of the columnar state is confirmed by 
the simulation results for the parallel Yukawa site model. 
The snapshot in Fig. 2a displays a hexagonal pattern in- 
dicative of freezing perpendicular to the nematic director. 
The columnar nature of the structure is refiected explic- 



itly in the pair correlation functions shown in Fig. 3. As 
to gi, the height of the first peak and the double-peaked 
shape of the second are a hallmark of long-ranged posi- 
tional order. More importantly, the flatness of gp signals 
an absence of long-range structure along the director (the 
small peak is due to weak mfra-columnar correlations). 
This means that there is no evidence of additional lay- 
ering in the z-direction indicative of three-dimensional 
crystalline order. To roughly estimate whether the state 
point adopted in the simulations corresponds to the bi- 
furcation diagram in Fig. 1 we may take the hydrody- 
namic aspect ratio Xh = 16.7 and thickness to estimate 
(kD)-^ ~ 2.4 and 77 - (7r/4)px^ = 0.011. Although the 
screening length is beyond the scale depicted in Fig. 1 
it can be easily inferred by extrapolation that the state 
point falls roughly in the columnar stability regime. 

Let us now turn to the case of freely rotating rods. The 
main question is whether or not the scenario sketched 
above is altered qualitatively if rotational degrees of free- 
dom are allowed for. The corresponding bifurcation dia- 
gram in Fig. 4 shows that this is not the case. In fact, the 
pre-emption of the nematic-columnar seems more out- 
spoken here: the crossover point has shifted to lower 
values and the the difference between the N-S and N- 
C volume fraction is enhanced beyond the intersection 
point. However, in order to be able to make a final as- 
sessment of the stability of the inhomogeneous phases 
we have to verify the stability of the nematic state it- 
self. To illustrate this, the bifurcation from the isotropic 
to the nematic state has been included in Fig. 4. This 
curve is obtained from a simple analysis described in Ap- 
pendix A. Below the I-N bifurcation, the nematic phase 
is no longer stable with respect to the isotropic and all 
instabilities of the nematic state located in this region 
become meaningless. For this reason, the scope of the 
bifurcation diagram is bounded by the intersection point 
at (kD)^^ sa 0.4. Below this point, the overall scenario is 
a nematic-smectic transition occurring at high screening 
{kD)~^ < 0.27 while a gradual crossover towards colum- 
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1 (K D)-^ 



FIG. 4: Bifurcation diagram for freely rotating charged rods. 
Plotted is the normalized Debye length kT^/D versus the 
hardcore volume fraction r]. The isotropic-nematic (I-N) bi- 
furcation is obtained from Eq. (|56|) in Appendix A. 



FIG. 5: Behavior of the Gaussian variational parameter a [see 
Eq. (go])] at the N-S (solid) and N-C (dotted) bifurcations as 
a function of Debye length. The inset shows the normalized 
a = a/oiH at the N-C bifurcation. 



nar order is expected at low salt conditions. The regime 
of low screening [i.e. {kD)~^ > 0.4], cannot be accessed 
but it is reasonable to anticipate translational symmetry- 
breaking instabilities (of either columnar or crystalline 
signature) of the isotropic phase. These are beyond the 
scope of the present calculations. The Gaussian varia- 
tional parameter at the instability is recorded in Fig. 5 
and its large value (the minimum is at a w 50) supports 
the use of the asymptotic analysis adopted throughout 
this paper. 

If the isotropic state were to be suppressed by applying 
a strong aligning external field, the nematic phase will be 
stable irrespective of density and pronounced columnar 
order can be expected in a large portion of the phase dia- 
gram, as we see in Fig. 4 (neglecting the TN curve) . This 
is illustrated by the simulation results. Judging from the 
snapshots in Fig. 2b and 2c it is clear that the columnar 
structures are robust against small orientational fluctua- 
tions. The differences between the layer pair correlation 
functions in Fig. 3 are very small. The number of lat- 
tice defects seems to increase slightly upon lowering the 
field strength but only at the smallest value (^ = 10) 
do we observe that the hexagonal pattern has vanished 
completely. The system nevertheless displays significant 
liquid like order perpendicular to the director which is 
typical for a dense nematic state. In the absence of the 
field, the systems is completely isotropic, in qualitative 
accordance with the phase diagram in Fig. 4. 

So far, we have implicitly assumed that the bifurca- 
tions represent thermodynamically stable phase transi- 
tions. Although the simulations clearly point to a (me- 
chanically) stable columnar phase for both parallel and 
asymptotically rotating rods, it would be desirable to get 
similar confirmation from the theory. To that end we 
have conducted a parametric expansion of the free en- 
ergy around the bifurcation, elaborated in Appendix B. 
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FIG. 6: Landau coefficient C4 given by Eq. (f65|l at the N-S 
(solid) and N-C (dotted) bifurcations as a function of Debye 
length corresponding to Fig. 4. 



The resulting Landau coefficient is depicted in Fig. 6. 
and it is clear that the thermodynamic stability condi- 
tion C4 < is generally fulfilled (except in the smectic 
region where N-S is pre-empted by N-C). A similar out- 
come is found for the parallel case, not shown here. 



IX. CONCLUSIONS 

An extensive bifurcation analysis is presented on trans- 
lational symmetry-breaking instabilities in nematic sys- 
tems of charged rods within Onsager's second virial 
theory. Starting from an artificial system of paral- 
lel charged rods, modelled via a hard-core Yukawa site 
model, nematic-smectic and nematic- (hexagonal) colum- 
nar instabilities are scrutinized as a function of the Debye 
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length. The treatment is then extended towards the more 
realistic situation of rods with rotational degrees of free- 
dom by employing an asymptotic Gaussian analysis of 
the orientation-dependent quantities. In both cases, dis- 
tinct preferential columnar order is observed for moder- 
ate Debye screening lengths. As a supplement, Brownian 
dynamics simulations for the point Yukawa site model 
were carried out and the stability of columnar struc- 
tures at low screening conditions and low particle con- 
centrations is qualitatively reproduced. No evidence for 
crystalline order was found. Based on this, no attempt 
has been made in theory to seek possible nematic-crystal 
instabilities, which could in principle be scrutinized by 
imposing coupled smectic and columnar density modu- 
lations. Moreover, these type of instabilities are more 
likely to occur at very low screening conditions where 
the effective rod aspect-ratio (which incorporates the ex- 
tent of the double layers) is no longer large enough to 
guarantee a stable nematic phase. In this region one 
may expect an isotropic fluid of rods to directly freeze 
into a crystalline lattice without intervention of nematic 
order^''^. These type of bifurcations are not addressed in 
this paper because the numerical effort involved in quan- 
tifying the necessary electrostatic end cap contributions 
for freely rotating rods is beyond the scale of the present 
calculations. 

Looking at the experimental results for charged rods 
we may put forward that a crossover from smectic or- 
der to a more intricate ordered state has been observed 
in concentrated systems of TMV rods upon decreasing 
ionic strengtbi-. However, at present it is not fully clear 
whether these structures are really columnar or represent 
three-dimensional crystalline order. More detailed struc- 
tural investigations are probably required to resolve this 
ambiguity. 

Our calculations also show that the columnar order 
can be realized at fairly low concentrations if the rods 
are rendered in near-parallel configurations by an exter- 
nal field. This could be any field that couples primarily 
to the rod orientations, for example magnetic, electric 
or shear flow fields. Finally, we remark that our results 
are connected to observations in other systems of rod-like 
mesogens with soft interactions. Columnar phases occur 
in complex systems of stiff polyelectrolytes like DNA^i^ 
and may be induced by an external magnetic field in 
systems of dipolar colloids^ and lath-shaped, goethite 
colloids^. Although other complicating features such as 
length polydispersity, dipole-dipole or non-uniform site- 
site interactions are at play in these systems, it is intrigu- 
ing to see that similar columnar structures may occur in 
a relatively simple model system for soft rods considered 
here. Future work could be aimed at finding out whether 
a scenario such as in Fig. 1 is qualitatively reproducible 
for other soft potentials. It is also desirable to enlarge the 
scope of the simulations such that a wide area of concen- 
trations and rod potentials can be covered and better 
comparisons with theory and experiments can be made. 



Appendix A: Isotropic-nematic bifurcation 

The bifurcation from the isotropic to the nematic state 
can be analytically derived starting from the stationarity 
condition Eq. ^ which we may rewrite as 

In /(f7) = fi + p I /(l^')'^(O; n')d^' (52) 



where the constant /i arises from the normalization con- 
dition of the ODF. The Mayer kernel at zero wavevector 
is related to the second virial coefficient of two charged 
rods [cf. Eq. dSll)]: 



$(0; Vt') = -2L'^D^B{n, n')\ sm-f{n, n')\ (53) 

The contributions of 0{LD^) do not depend on orienta- 
tion and hence drop out of the free energy minimization 
with respect to the ODF leading to Eq. ([5^ . Close to 
the bifurcation point the ODF can be parametrized as: 



£7^2 (cos 0)] 



(54) 



with Vn a Legendre polynomial and e an arbitrarily small 
order parameter quantifying a uniaxial nematic pertur- 
bation to the isotropic ODF / = 1/An. Likewise, the 
orientation-dependent functions in the Mayer kernel can 
be expanded in terms of these polynomials in the follow- 



ing way^: 



F(sin7) = ca + C2V2{cose)V2{cose') + 



(55) 



with coefficients cq = it /A and C2 
and Co = (7r/4)[ln2- 1/2] and C2 
for F{x) = —x\nx. Inserting 
arity condition Eq. (f52|) . linear 
and some rearranging (using the 
^\d{cose)Vn[cose)V^{cose) = 2(5„,„/(2n+ l)) readily 
leads to a closed expression for the bifurcation concen- 
tration: 



= -57r/32 for F{x) = x 
= (-57r/32)[ln2-5/4] 
these into the station- 
izing with respect to e 
orthogonality condition 



L 



^ ^ D 1 + (Ki:))-i(lnA-t-7£; - ln2 



with A given by Eq. ([28 



(56) 



Appendix B: Landau expansion around the 
bifurcation point 

An expression for the free energy difference AF = 
Fi — Fm between the inhomogeneous state and the ho- 
mogeneous nematic reference state can be obtained by 
inserting the parametrization Eq. ^ into the functional 
Eq. (IT|). After some rearranging one can show that the 
free energy difference per particle reads: 



A/3f 
TV 



1 

2^ 



J dcmc)inw^(c)-f 



i^fik) (57) 
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using the short-hand notation ^f{lq) 
and 



wio 



cos(ZC) 



(58) 



1=1 



For the sake of brevity we will restrict ourselves to the 
smectic case here, i.e. a single instability mode. The 
analysis for the columnar symmetry can be carried out 
in a similar manner. Minimizing with respect to the or- 
der parameters ai and the wave number q leads to the 
following stability conditions 



P) 



1 

P 



2-iT 



dC cos(toC) In W{C,) 



am^f{mq) =0, m > 1 



pY,afmiq) = 

1=1 



(59) 



with $' — d<t/dq. Let us now propose the following 
expansions in terms of an arbitrarily small parameter e 



W{C;e) 
Pie) 



1 + eai cos C + £ cos 2^ -I- e 03 cos 3C - 
Po + epi 



2 

e P2 



q{e) = qo + eqi + ("q^ + 



(60) 



The expansion of W is justified close to bifurcation point 
where the spatial inhomogeneity is supposed to be weak. 
The zeroth order solution e = reproduces the bifur- 
cation condition po^f(q) = 1 and hence {ai, po, qo} = 
{1, p* , q*}. Inserting the parametrization into Eq. (j59|) 
allows us to expand the set of stationarity conditions as 
follows: 



A{e) = e(A(o)+eA«+e24(2,) + ...)=0 
B{e) = ( 



B 



(0) 



^2^(2) 



TO > 1 

= (61) 



The calculation of the coefficients is straightforward but 
tedious and we will only give the essential results in the 
remainder of the Appendix. Performing an order by or- 
der solution of the above set of equations yields up to 
first order 



Pi 







qi 

0.2 



4(i-po$/(2g)) 



(62) 



and up to second order (dropping the caret and the sub- 
script / for notational clarity) 



P2 



92 



l^"{q)il-2a2)+Poal<^'iqWi2q) 

i$(g)$"((Z) - i^'{q)r 
^<i>'{q){l-2a2) + al<i>'i2q) 

2^2 - 12 



= p„$(3,)-l . 
The free energy difference can now be expanded in an 
analogous way using the parametrization Eq. (j60[) in Eq. 
([57|) . The Landau free energy as a function of the order 
parameter e reads 



A(3F 
N 



e'C2 



(64) 



with coefficients 

C2 
C3 

C4 



/Oo^(2g) 
64 Vl-po$(2g) 



(65) 



In agreement with Mulder's results^^ we get a zero cu- 
bic contribution in Eq. (j64p indicating the transition 
towards the inhomogeneous state to be of second order. 
The thermodynamic stability of the new state however 
is guaranteed by the condition C4 < which has to be 
assessed numerically via the Mayer kernel. 
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